# install packages
install.packages("statnet")
install.packages("xergm")
install.packages("texreg")
install.packages("xtable")
install.packages("ggplot2")

# set directories
setwd("~/Dropbox/ISQ_ms/5 Final/Replication")
SAVEDIR <- "~/Dropbox/ISQ_ms/5 Final/Replication"

# load libraries
detach(package:igraph, unload=T)
require("statnet")
require("xergm")
require("texreg")
require("xtable")
library(ggplot2)

### Load data

# Read in edgelists
el1995 <- read.csv("mod4_el1995.csv",stringsAsFactors=F)
el2000 <- read.csv("mod4_el2000.csv",stringsAsFactors=F)
el2005 <- read.csv("mod4_el2005.csv",stringsAsFactors=F)

# Read in dyadic covariates

# alliance networks
ally1995 <- read.csv("mod4_ally1995.csv",row.names=1)
ally1995 <- as.matrix(ally1995)
ally2000 <- read.csv("mod4_ally2000.csv",row.names=1)
ally2000 <- as.matrix(ally2000)
ally2005 <- read.csv("mod4_ally2003.csv",row.names=1)
ally2005 <- as.matrix(ally2005)

# contiguity matrix
cont1995 <- read.csv("mod4_cont1995.csv",row.names=1)
cont1995 <- as.matrix(cont1995)
cont2000 <- read.csv("mod4_cont2000.csv",row.names=1)
cont2000 <- as.matrix(cont2000)
cont2005 <- read.csv("mod4_cont2005.csv",row.names=1)
cont2005 <- as.matrix(cont2005)

# trade network
trade1995 <- read.csv("mod4_trade1995.csv",row.names=1)
trade1995 <- as.matrix(trade1995)
trade2000 <- read.csv("mod4_trade2000.csv",row.names=1)
trade2000 <- as.matrix(trade2000)
trade2005 <- read.csv("mod4_trade2005.csv",row.names=1)
trade2005 <- as.matrix(trade2005)

# Read in vertex covariates
vcov <- read.csv("mod4_vcov.csv",stringsAsFactors=F)
vcov1995 <- vcov[vcov$year==1995,]
vcov2000 <- vcov[vcov$year==2000,]
vcov2005 <- vcov[vcov$year==2005,]


### Create network

# Initialize networks and assign labels

net1995 <- network.initialize(nrow(vcov1995))
network.vertex.names(net1995) <- vcov1995$stateabb

net2000 <- network.initialize(nrow(vcov2000))
network.vertex.names(net2000) <- vcov2000$stateabb

net2005 <- network.initialize(nrow(vcov2005))
network.vertex.names(net2005) <- vcov2005$stateabb

# Add in edges
net1995[as.matrix(el1995)] <- 1
net2000[as.matrix(el2000)] <- 1
net2005[as.matrix(el2005)] <- 1

# Define vertex attributes

# For 1995

set.vertex.attribute(net1995,"pop", vcov1995$log_pop)
set.vertex.attribute(net1995,"gdp", vcov1995$log_realgdp)
set.vertex.attribute(net1995,"gdppc", vcov1995$log_rgdppc)

set.vertex.attribute(net1995,"cinc", vcov1995$log_cinc)
set.vertex.attribute(net1995,"milex", vcov1995$log_milex)
set.vertex.attribute(net1995,"milper", vcov1995$log_milper)

set.vertex.attribute(net1995,"nw", vcov1995$nw)

set.vertex.attribute(net1995,"polity", vcov1995$polity)
set.vertex.attribute(net1995,"dem_dum", vcov1995$dem_dum)

set.vertex.attribute(net1995,"pts_inv", vcov1995$pts_inv)
set.vertex.attribute(net1995,"pts_dum", vcov1995$pts_dum)

set.vertex.attribute(net1995,"hf_efiscore", vcov1995$hf_efiscore)
set.vertex.attribute(net1995,"ief_dum", vcov1995$ief_dum)

set.vertex.attribute(net1995,"region", vcov1995$sub_region_n)
set.vertex.attribute(net1995,"igo_hq", vcov1995$igo_hq)

# For 2000

set.vertex.attribute(net2000,"pop", vcov2000$log_pop)
set.vertex.attribute(net2000,"gdp", vcov2000$log_realgdp)
set.vertex.attribute(net2000,"gdppc", vcov2000$log_rgdppc)

set.vertex.attribute(net2000,"cinc", vcov2000$log_cinc)
set.vertex.attribute(net2000,"milex", vcov2000$log_milex)
set.vertex.attribute(net2000,"milper", vcov2000$log_milper)

set.vertex.attribute(net2000,"nw", vcov2000$nw)

set.vertex.attribute(net2000,"polity", vcov2000$polity)
set.vertex.attribute(net2000,"dem_dum", vcov2000$dem_dum)

set.vertex.attribute(net2000,"pts_inv", vcov2000$pts_inv)
set.vertex.attribute(net2000,"pts_dum", vcov2000$pts_dum)

set.vertex.attribute(net2000,"hf_efiscore", vcov2000$hf_efiscore)
set.vertex.attribute(net2000,"ief_dum", vcov2000$ief_dum)

set.vertex.attribute(net2000,"region", vcov2000$sub_region_n)
set.vertex.attribute(net2000,"igo_hq", vcov2000$igo_hq)

# For 2005

set.vertex.attribute(net2005,"pop", vcov2005$log_pop)
set.vertex.attribute(net2005,"gdp", vcov2005$log_realgdp)
set.vertex.attribute(net2005,"gdppc", vcov2005$log_rgdppc)

set.vertex.attribute(net2005,"cinc", vcov2005$log_cinc)
set.vertex.attribute(net2005,"milex", vcov2005$log_milex)
set.vertex.attribute(net2005,"milper", vcov2005$log_milper)

set.vertex.attribute(net2005,"nw", vcov2005$nw)

set.vertex.attribute(net2005,"polity", vcov2005$polity)
set.vertex.attribute(net2005,"dem_dum", vcov2005$dem_dum)

set.vertex.attribute(net2005,"pts_inv", vcov2005$pts_inv)
set.vertex.attribute(net2005,"pts_dum", vcov2005$pts_dum)

set.vertex.attribute(net2005,"hf_efiscore", vcov2005$hf_efiscore)
set.vertex.attribute(net2005,"ief_dum", vcov2005$ief_dum)

set.vertex.attribute(net2005,"region", vcov2005$sub_region_n)
set.vertex.attribute(net2005,"igo_hq", vcov2005$igo_hq)


# Check networks
net1995
net2000
net2005


### For the 1995-2005 period

# t=1 -> 1995
# t=2 -> 2000
# t=3 -> 2005


# Combine dependent networks 

all_dip_nets <- list(as.matrix(net1995), as.matrix(net2000), as.matrix(net2005))

# Combine dyadic covariates

all_ally_nets <- list(ally1995, ally2000)
all_cont_nets <- list(cont1995, cont2000)
all_trade_nets <- list(trade1995, trade2000)


# check objects
class(all_dip_nets)
length(all_dip_nets)

# Use the preprocess function on the outcome network
# Tell it what missing values and structural zeros (instances in which tie impossible) look like

# estimation (dep nets) starts at t=2 and ends at t=3
all_dep <- preprocess(all_dip_nets,
                      lag = TRUE, covariate = FALSE, na = NA,
                      na.method = "fillmode", structzero = 10,
                      structzero.method = "remove") 

# Create lagged dependent network
all_mem <- preprocess(all_dip_nets, 
                      lag = TRUE, covariate = TRUE, memory = "autoregression",
                      na = NA, na.method = "fillmode", structzero = 10,
                      structzero.method = "remove")

# Check that everything looks right
length(all_dep)
sapply(all_dep, dim)
sapply(all_mem, dim)
# dim of the mem term need to be the same (preprocess does that)

# Preprocess dyadic covariates
# covariates start at t= 1 and end at t=5

all_ally <- preprocess(all_ally_nets, all_dep, 
                       lag = FALSE, covariate = TRUE)
str(all_ally)
sapply(all_ally, dim)

all_cont <- preprocess(all_cont_nets, all_dep, 
                       lag = FALSE, covariate = TRUE)
str(all_cont)
sapply(all_cont, dim)

all_trade <- preprocess(all_trade_nets, all_dep, 
                        lag = FALSE, covariate = TRUE)
str(all_trade)
sapply(all_trade, dim)


## Set vertex attributes

# get vertex names for the dependent networks
str(all_dep[[1]])
stateIDs2000 <- colnames(all_dep[[1]])
stateIDs2005 <- colnames(all_dep[[2]])

# turn dep into networks to be able to set vertex attributes
for (i in 1:length(all_dep)) {
  all_dep[[i]] <- network(all_dep[[i]])}
all_dep[[1]]


# For 1995 covariates (2000 dependent net)

pop1995 <- vcov1995$log_pop[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"pop", pop1995)

gdp1995 <- vcov1995$log_realgdp[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"gdp", gdp1995)

gdppc1995 <- vcov1995$log_rgdppc[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"gdppc", gdppc1995)

log_milex1995 <- vcov1995$log_milex[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"milex", log_milex1995)

log_milper1995 <- vcov1995$log_milper[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"milper", log_milper1995)

log_cinc1995 <- vcov1995$log_cinc[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"cinc", log_cinc1995)

nw1995 <- vcov1995$nw[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"nw", nw1995)

polity1995 <- vcov1995$polity[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"polity", polity1995)

dem_dum1995 <- vcov1995$dem_dum[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"dem_dum", dem_dum1995)

pts_inv1995 <- vcov1995$pts_inv[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"pts_inv", pts_inv1995)

pts_dum1995 <- vcov1995$pts_dum[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"pts_dum", pts_dum1995)

hf_efiscore1995 <- vcov1995$hf_efiscore[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"hf_efiscore", hf_efiscore1995)

ief_dum1995 <- vcov1995$ief_dum[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"ief_dum", ief_dum1995)

region1995 <- vcov1995$sub_region_n[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"region", region1995)

igo_hq1995 <- vcov1995$igo_hq[vcov1995$stateabb %in% stateIDs2000]
set.vertex.attribute(all_dep[[1]],"igo_hq", igo_hq1995)

# check
all_dep[[1]]


# For 2000 covariates (2005 dependent net)

pop2000 <- vcov2000$log_pop[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"pop", pop2000)

gdp2000 <- vcov2000$log_realgdp[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"gdp", gdp2000)

gdppc2000 <- vcov2000$log_rgdppc[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"gdppc", gdppc2000)

log_milex2000 <- vcov2000$log_milex[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"milex", log_milex2000)

log_milper2000 <- vcov2000$log_milper[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"milper", log_milper2000)

log_cinc2000 <- vcov2000$log_cinc[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"cinc", log_cinc2000)

nw2000 <- vcov2000$nw[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"nw", nw2000)

polity2000 <- vcov2000$polity[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"polity", polity2000)

dem_dum2000 <- vcov2000$dem_dum[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"dem_dum", dem_dum2000)

pts_inv2000 <- vcov2000$pts_inv[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"pts_inv", pts_inv2000)

pts_dum2000 <- vcov2000$pts_dum[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"pts_dum", pts_dum2000)

hf_efiscore2000 <- vcov2000$hf_efiscore[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"hf_efiscore", hf_efiscore2000)

ief_dum2000 <- vcov2000$ief_dum[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"ief_dum", ief_dum2000)

region2000 <- vcov2000$sub_region_n[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"region", region2000)

igo_hq2000 <- vcov2000$igo_hq[vcov2000$stateabb %in% stateIDs2005]
set.vertex.attribute(all_dep[[2]],"igo_hq", igo_hq2000)

# check
all_dep[[2]]


###################################### Estimate models ##############################################

set.seed(10) 
# Main model (Model 1)
model4e <- btergm(all_dep ~ edges + istar(2) + ostar(2) + mutual + triangles
                  +absdiff("polity")
                  +absdiff("pts_inv")
                  +absdiff("hf_efiscore")
                  +absdiff("gdppc")
                  +absdiff("milex")
                  +absdiff("nw")
                  +nodeicov("polity")
                  +nodeicov("pts_inv")
                  +nodeicov("hf_efiscore")
                  +nodeicov("gdppc")
                  +nodeicov("milex")
                  +nodeicov("nw")
                  +edgecov(all_ally)
                  +edgecov(all_trade)
                  +edgecov(all_cont)
                  +nodematch("region")
                  +nodeicov("igo_hq")
                  +nodeocov("gdppc")
                  +edgecov(all_mem),
                  R=1000)


### Figure 3. Temporal Exponential Random Graph Model of Diplomatic Ties for 1995-2005

# Put model estimates into temporary data.frames:
model4eFrame <- data.frame(Variable = c("edges", "Popularity", "Sociality", "Reciprocity", "Transitivity",
                                        "Difference in Democracy",
                                        "Difference in Human Rights",
                                        "Difference in Economic Freedom",
                                        "Difference in GDP/Capita", "Difference in Military Spending", 
                                        "Difference in Nuclear Weapons", 
                                        "Democracy(j)", "Human Rights(j)", "Economic Freedom(j)", 
                                        "GDP/Capita(j)", "Military Expenditure(j)", "Nuclear Weapons(j)",
                                        "Alliance", "Trade", "Contiguity", "Same Region",
                                        "IGO Headquarters", "GDP/Capita(i)", "Tie Stability"),
                           Coefficient = summary(model4e)[,1],
                           CIL = summary(model4e)[,2],
                           CIU = summary(model4e)[,3])

# make var order the same as in the model
model4eFrame <- model4eFrame[-1,] # delete edges coeff
model4eFrame$Variable <- factor(model4eFrame$Variable, levels = unique(rev(model4eFrame$Variable)))

# create sig variable to give different colors to significant/insignificant coeffs
model4eFrame$sig <- 1
model4eFrame$sig[model4eFrame$CIL<0 & model4eFrame$CIU>0] <- 0
model4eFrame$sig <- factor(model4eFrame$sig)

# Plot (B&W)
pdf(file=file.path(SAVEDIR, "DuqueFig3.pdf"), width=6.5, height=5.2)
p1a <- ggplot(model4eFrame, aes(x=Variable, fill = sig, colour = sig)) 
p1a <- p1a + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
p1a <- p1a + geom_linerange(aes(ymin = CIL, ymax = CIU), colour = "black")
p1a <- p1a + geom_point(aes(x = Variable, y = Coefficient), shape = 21, size = 1.2)
p1a <- p1a + coord_flip() + theme_bw()
p1a <- p1a + scale_fill_manual(values = c("white", "black"))
p1a <- p1a + scale_colour_manual(values = c("black", "black"))
p1a <- p1a + geom_vline(xintercept=c(7.5, 13.5, 19.5), size=.25)
p1a <- p1a + theme(legend.position="none")
p1a <- p1a + annotate("text", x = 2.5, y = 3, label = "Controls", size=2.7)
p1a <- p1a + annotate("text", x = 8.5, y = 3, label = "State attributes", size=2.7)
p1a <- p1a + annotate("text", x = 14.5, y = 3, label = "Homophily", size=2.7)
p1a <- p1a + annotate("text", x = 20.5, y = 3, label = "Endogenous effects", size=2.7)
p1a <- p1a + theme(axis.text.y = element_text(size=7.6))
p1a <- p1a + theme(axis.text.x = element_text(size=7.6))
p1a <- p1a + scale_y_continuous(name="")
p1a <- p1a + scale_x_discrete(name="", breaks=levels(model4eFrame$Variable),
                              labels=c("Tie Stability",expression(GDP~capita[i]),expression(IGO~Headquarters[j]),
                                       "Same Region","Contiguity",
                                       "Trade","Alliance",
                                       expression(Nuclear~Weapons[j]),expression(Military~Spending[j]),
                                       expression(GDP~per~capita[j]),expression(Economic~Freedom[j]),
                                       expression(Human~Rights[j]),expression(Democracy[j]),
                                       expression(Nuclear~Weapons[list(abs(i-j))]),
                                       expression(Military~Spending[list(abs(i-j))]), 
                                       expression(GDP~per~capita[list(abs(i-j))]),
                                       expression(Economic~Freedom[list(abs(i-j))]),
                                       expression(Human~Rights[list(abs(i-j))]),
                                       expression(Democracy[list(abs(i-j))]),
                                       "Transitivity","Reciprocity","Sociality","Popularity"))
p1a
dev.off()


# Substantive interpretation (for statistically significant coefficients only)
# endogenous effects 
exp(summary(model4e)[[2]]*(51-23)) # pop (comparison btw Cuba and Bolivia)
exp(summary(model4e)[[4]]) # recip
exp(summary(model4e)[[5]]*(1642-1615)) # transitivity (comparison btw HUN->LIB and HUN->CHA)


### Table 3. Change in Odds of Tie from Typical Change in Variable

change <- sd(abs(apply(combn(polity2000,2), 2, diff))) #Std. deviation for absdiff polity
exp(coef(model4e)[['absdiff.polity']]*(-change))

change <- IQR(abs(apply(combn(pts_inv2000,2), 2, diff))) #IQR for absdiff pts_inv (categorical)
exp(coef(model4e)[['absdiff.pts_inv']]*(-change))

change <- sd(abs(apply(combn(hf_efiscore2000,2), 2, diff)))
exp(coef(model4e)[['absdiff.hf_efiscore']]*(-change))

change <- sd(abs(apply(combn(gdppc2000,2), 2, diff)))
exp(coef(model4e)[['absdiff.gdppc']]*(-change))

change <- sd(abs(apply(combn(log_milex2000,2), 2, diff)))
exp(coef(model4e)[['absdiff.milex']]*(-change))

exp(coef(model4e)[['nodeicov.polity']]*(sd(polity2000)))
exp(coef(model4e)[['nodeicov.pts_inv']]*(IQR(pts_inv2000)))
exp(coef(model4e)[['nodeicov.hf_efiscore']]*(sd(hf_efiscore2000)))
exp(coef(model4e)[['nodeicov.gdppc']]*(sd(gdppc2000)))
exp(coef(model4e)[['nodeicov.milex']]*(sd(log_milex2000)))
exp(coef(model4e)[['nodeicov.nw']]*(1))



### Figure A3. Goodness-of-Fit of the TERGM of Diplomatic Ties for 1995-2005 (Model 1)

# goodness of fit
set.seed(10)
gof4e <- gof(model4e, nsim = 1000) 
plot(gof4e) 

# save GOF
pdf(file=file.path(SAVEDIR, "DuqueFigA3.pdf"), width=10, height=8)
plot(gof4e, boxplot.mfrow=T, roc=F, pr=F)
dev.off()



### Table A12. Degeneracy Checks for the TERGM of Diplomatic Ties for 1995-2005 (Model 1)

# degeneracy check
set.seed(10)
mod4edegen <- checkdegeneracy(model4e, nsim = 1000)
degentab <- cbind(as.data.frame(mod4edegen[1]), as.data.frame(mod4edegen[2]))

# save
require("xtable")
x.long <- xtable(degentab,
                 caption = "Degeneracy Checks for the TERGM of Diplomatic Ties for 1995-2005",
                 label = "Degeneracy1995-2005",
                 align = "lrrrrrrrrrrrr")
rownames(x.long) <- c("Edges", "Popularity", "Sociality", "Reciprocity", "Transitivity",
                      "Democracy$_{|i-j|}$","Human Rights$_{|i-j|}$","Economic Freedom$_{|i-j|}$",
                      "GDP/capita$_{|i-j|}$", "Military Spending$_{|i-j|}$", "Nuclear Weapons$_{|i-j|}$", 
                      "Democracy$_j$", "Human Rights$_j$", "Economic Freedom$_j$", 
                      "GDP/Capita$_j$", "Military Expenditure$_j$", "Nuclear Weapons$_j$",
                      "Alliance", "Trade", "Contiguity", "Same Region",
                      "IGO Headquarters", "GDP/capita$_i$", "Tie Stability")
print(x.long, floating.environment = 'sidewaystable',
      include.rownames = T)



### Table A13. Temporal Exponential Random Graph Models of Diplomatic Ties for 1995-2005 (Robustness Checks)

set.seed(10)
# Model 2: logit with node attributes only (no homophily)
model4b <- btergm(all_dep ~ edges 
                  +nodeicov("polity")
                  +nodeicov("pts_inv")
                  +nodeicov("hf_efiscore")
                  +nodeicov("gdppc")
                  +nodeicov("milex")
                  +nodeicov("nw")
                  +edgecov(all_ally)
                  +edgecov(all_trade)
                  +edgecov(all_cont)
                  +nodematch("region")
                  +nodeicov("igo_hq")
                  +nodeocov("gdppc")
                  +edgecov(all_mem),
                  R=1000)

set.seed(10)
# Model 3: logit with homophily
model4c <- btergm(all_dep ~ edges
                  +absdiff("polity")
                  +absdiff("pts_inv")
                  +absdiff("hf_efiscore")
                  +absdiff("gdppc")
                  +absdiff("milex")
                  +absdiff("nw")
                  +nodeicov("polity")
                  +nodeicov("pts_inv")
                  +nodeicov("hf_efiscore")
                  +nodeicov("gdppc")
                  +nodeicov("milex")
                  +nodeicov("nw")
                  +edgecov(all_ally)
                  +edgecov(all_trade)
                  +edgecov(all_cont)
                  +nodematch("region")
                  +nodeicov("igo_hq")
                  +nodeocov("gdppc")
                  +edgecov(all_mem),
                  R=1000)

texreg(c(model4b, model4c, model4e),
       digits=3, single.row=T, stars = 0.05, sideways=T)



### Table A15. Temporal Exponential Random Graph Models of Diplomatic Ties for 1995-2005 (Robustness Checks)

set.seed(10)
# Model 8: replaces per capita GDP with GDP & pop
model4d <- btergm(all_dep ~ edges + istar(2) + ostar(2) + mutual + triangles
                  +absdiff("polity")
                  +absdiff("pts_inv")
                  +absdiff("hf_efiscore")
                  +absdiff("gdp")
                  +absdiff("pop")
                  +absdiff("milex")
                  +absdiff("nw")
                  +nodeicov("polity")
                  +nodeicov("pts_inv")
                  +nodeicov("hf_efiscore")
                  +nodeicov("gdp")
                  +nodeicov("pop")
                  +nodeicov("milex")
                  +nodeicov("nw")
                  +edgecov(all_ally)
                  +edgecov(all_trade)
                  +edgecov(all_cont)
                  +nodematch("region")
                  +nodeicov("igo_hq")
                  +nodeocov("gdp")
                  +edgecov(all_mem),
                  R=1000)

set.seed(10)
# Model 9: replaces milex with CINC score 
model4z <- btergm(all_dep ~ edges + istar(2) + ostar(2) + mutual + triangles
                  +absdiff("polity")
                  +absdiff("pts_inv")
                  +absdiff("hf_efiscore")
                  +absdiff("gdppc")
                  +absdiff("cinc")
                  +absdiff("nw")
                  +nodeicov("polity")
                  +nodeicov("pts_inv")
                  +nodeicov("hf_efiscore")
                  +nodeicov("gdppc")
                  +nodeicov("cinc")
                  +nodeicov("nw")
                  +edgecov(all_ally)
                  +edgecov(all_trade)
                  +edgecov(all_cont)
                  +nodematch("region")
                  +nodeicov("igo_hq")
                  +nodeocov("gdppc")
                  +edgecov(all_mem),
                  R=1000)

texreg(c(model4d, model4z, model4e), digits=3, single.row=T, sideways=T)

